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We present a numerical study of magnetic ordering in spin ice on kagome, a two- 
^ I dimensional lattice of corner-sharing triangles. The magnet has six ground states 

and the ordering occurs in two stages, as one might expect for a six-state clock 
model. In spin ice with short-range interactions up to second neighbors, there is 
an intermediate critical phase separated from the paramagnetic and ordered phases 
by Kosterlitz-Thouless transitions. In dipolar spin ice, the intermediate phase has 
long-range order of staggered magnetic charges. The high and low-temperature 
phase transitions are of the Ising and 3-state Potts universality classes, respec- 
tively. Freeze-out of defects in the charge order produces a very large spin correla- 
tion length in the intermediate phase. As a result of that, the lower-temperature 
transition appears to be of the Kosterlitz-Thouless type. 
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1. Introduction 

Spin ice is a ferromagnet with peculiar properties (Gingras 2011). Its discovery in 
the rare-earth pyrochlore IIo2Ti207 by Harris et al. (1997) attracted interest of 
physicists because geometric frustration results in an enormous number of ground 
states in spin ice. The degeneracy is similar to that of proton positions in water 
ice (Petrenko & Whitworth 1999) and manifests itself in a large residual entropy 
at low temperatures measured by Ramirez et al. (1999). Later it was shown that 
spin ice has another peculiar feature: even though it remains disordered as the 
temperature goes to zero, spin correlations decay with the distance algebraically 
rather than exponentially (Isakov et al. 2004; Henley 2005). Such critical behavior 
is characteristic of systems with constraints exemplified by lattice models with 
hardcore dimers. In the case of spin ice, at low temperatures each tetrahedron 
of rare-earth ions is forced to have two spins pointing toward the center of the 
tetrahedron and two away from it. 

The current resurgence of interest in spin ice stems from an unusual charac- 
ter of magnetic excitations in this class of materials (Ryzhkin 2005; Castelnovo, 
Moessner & Sondhi 2008). Reversing a single spin in a spin- ice state violates the 
above-mentioned constraint on the two tetrahedra sharing that spin. One of them 
has three spins pointing in and one pointing out, the other one spin pointing in and 
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Figure 1. (a) Generic spin-ice state on kagome. (b) A spin-ice (micro)state with magnetic 
charges ordered in a staggered patter, (c) Expectation values of spins and magnetic charges 
in a charged-ordered intermediate state. Note that the expectation values for the Ising 
variables are (at) = ±1/3 in this state, so the charges on triangles are (Qa) ~ ±1. (d) A 
ground state of kagome spin ice with ordered charges and magnetic dipoles. The dashed 
line is the perimeter of a magnetic unit cell. 

three out. By reversing additional spins, the two defects can be moved around in- 
dependently from each other. These two defect tetrahedra carry equal and opposite 
magnetic charges defined as sources and sinks of magnetic field H. A defect with 
magnetic charge Q in an external magnetic field Hoxt experiences a Zeeman force 
/ioQHoxt- Two defects with magnetic charges Qi and Q2 interact with one another 
via a Coulomb potential /ioQiQ2/(47rr). Various properties of spin ice can be most 
naturally described by focusing on the motion of these defects (Fennell et al. 2009; 
Jaubert et al. 2009; Morris et al. 2009; Ryzhkin & Ryzhkin 2011). 

The concept of magnetic charge is not exactly novel. It can be found in earlier 
research articles (Saitoh et al. 2004) and even in textbooks (Landau & Lifshitz 
1984). Magnetic charges in spin ice are remarkable because they are mobile and 
represent a rare example of fractionalized excitations in three spatial dimensions: 
the underlying degrees of freedom, magnetic dipoles, must be split in half, so to 
speak, to create magnetic monopoles. 

It is worth pointing out that magnetic charges in spin ice are fundamentally 
different from the magnetic monopoles introduced by Dirac (1931) to explain quan- 
tization of electric charge. They are sources and sinks of magnetic field H, rather 
than of magnetic induction B. For that reason, there are no restrictions on the 
possible values of magnetic charges in spin ice. 

In this paper, we illustrate the utility of the concept of magnetic charge for the 
problem of spin ice on a different lattice. A natural generalization of the pyrochlore 
network is a lattice consisting of corner-sharing simplexes (Chalker 2011). In two 
dimensions, corner-sharing triangles form kagome. Fig. [TJ A spin is shared by two 
triangles and points along the line connecting their centers. The easy axes of the 
three spins on a triangle make angles of 120° with each other. Like on the pyrochlore 
lattice, interactions of the antiferromagnetic kind are not frustrating: an isolated 
triangle has two ground states with all spins pointing in or all pointing out. A 
kagome lattice with nearest-neighbor antiferromagnetic interactions also has two 
ground states, in which spins point into triangles of the same orientation. On the 
other hand, ferromagnetic interactions are frustrating and yield six ground states 
on a single triangle, with two spins pointing in (and one out) or two pointing out 
(and one in). The number of ground states grows exponentially with the number 
of spins when interactions are ferromagnetic and only include nearest neighbors. 
In that, spin ice on kagome is quite similar to its analog on the pyrochlore lattice. 
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However, the addition of longer-range interactions, especially dipolar ones, reveals 
some important differences between the two systems. 

The main difference between the pyrochlorc network and kagome concerns the 
allowed values of magnetic charge Q of a simplex defined as the flux of magnetiza- 
tion M into the simplex. For a tetrahedron, Q is even (in the appropriate units): 
zero in a two-in-two-out state, ±2 in thrcc-in-one-out or three-out-onc-in states and 
±4 in the all-in or all-out states. The energy of spin ice penalizes states with a large 
(in absolute value) charge Q so that simplexes are neutral, Q = 0, in low-energy 
states. The lack of magnetic charges on tetrahedra explains a piizzling observation 
by Gingras and coworkers: spin-ice states remain very nearly degenerate even in 
the presence of long-range dipolar interactions (Gingras & den Hertog 2001). As 
Castelnovo et al. (2008) showed, the energy of dipolar interactions is well approxi- 
mated by the Coulomb energy of magnetic charges residing on simplexes. Because 
these charges vanish in two-in-two-out states, their Coulomb energy is the same. 
The residual interactions, neglected in the Coulomb approximation, force the sys- 
tem into a phase with long-range magnetic order at temperatures low compared to 
the strength of dipolar interactions (Melko et al. 2001). 

In contrast, the magnetic charge of a triangle is ±1 in a spin-ice state and ±3 
when its spins point all in or all out. Even if the interactions tend to suppress 
magnetic charges, the lowest (in magnitude) values of charge on a triangle are ±1. 
This has important consequences that we briefly outline below and consider in 
detail in the rest of the paper. 

First, the presence of a charge degree of freedom on simplexes makes spin ice 
on kagome much more fluid than its pyrochlorc counterpart. In pyrochlorc spin 
ice, spin dynamics slows down considerably at low temperatures. Single-spin flips 
result in the creation of two magnetic charges and thus become forbidden when the 
temperature falls below the energy cost of a monopole. Moves within the low-energy 
sector require flipping entire chains of alternating spins with a minimal length of 
six (Melko et al. 2001). Alternatively, spin fluctuations can be mediated by the 
motion of a few remaining magnetic monopolcs. Experimental studies confirm a 
very slow relaxation in spin ice at liquid- helium temperatures (Snyder et al. 2001; 
Giblin et al. 2011), although the primary mechanism of low-T spin dynamics in 
the rare-earth pyrochlores remains to be clarified. No such impediment to spin 
fluctuations exists in kagome ice because a single spin flip does not necessarily 
take the magnet out of the low-energy sector. Such a local move is allowed if the 
magnetic charges of the two simplexes sharing the spin change from +1 and —1 to 
— 1 and -|-1. Although at present there are no materials realizing two-dimensional 
kagome spin ice, several experimental groups have made artiflcial magnetic arrays 
with this geometry (Tanaka et al. 2006). Qi et al. (2008) demonstrated that their 
artificial spin ice on kagome stayed strictly within the low-energy sector in which 
the magnetic charges remain ±1. Ladak et al. (2010) observed triple charges, but 
this is likely the result of strong quenched disorder in their samples (Ladak et al. 
2011). 

Second, the presence of a charge degree of freedom opens a possibility of new 
kinds of ordered phases in spin ice. Pyrochlorc spin ice is expected to have only two 
thermodynamics phases. As the temperature is lowered from -l-oo, the paramagnetic 
state gradually crosses over to the spin-ice state, which still preserves all of the 
symmetries of the Hamiltonian and thus remains a paramagnet, albeit a correlated 
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one, from the standpoint of the Landau theory. At a sufficiently low temperature, 
the residual spin interactions not captured by the Coulomb model (Castelnovo et 
al. 2008) induce a phase transition into an ordered state that breaks the time- 
reversal and lattice symmetries (Melko & Gingras 2004). On kagome ice, one has 
good reasons to expect an intermediate ordered phase with staggered magnetic 
charges that minimizes the Coulomb energy (Mollcr and Moessncr 2009; Chcrn et 
al. 2011). This phase is characterized by very strong spin fluctuations. At the lowest 
temperatures, the symmetry is broken further and the system settles into a state 
with long-range magnetic order. 

To test these heuristic considerations, we have performed extensive numerical 
studies of spin ice on kagome, both with short-range interactions (nearest and next- 
nearest neighbors) and with long-range dipolar interactions. The problem turned 
out to be subtle and the answers different for the models with short and long-range 
interactions. In both cases, there is an intermediate phase. However, the nature 
of the intermediate phase is drastically different. The dipolar spin ice has the ex- 
pected charge-ordered phase, whereas the short-range model has a critical interme- 
diate phase without true long-range order. The determination of the universality 
classes of the phase transitions required rather large system sizes because some of 
the transitions were of the Kosterlitz-Thouless type. Simulating large systems is 
particularly difficult in the presence of long-range interactions. 

2. Kagome spin ice: the model 

(a) The Hamiltonians 

The model of spin ice on kagome was first introduced by Wills et al. (2002). 
It has Ising spins Sj = aiSi living on the vertices of corner-sharing triangles of a 
kagome lattice. The unit vector ej points along the line connecting the centers of 
the two triangles, while the Ising variable ai encodes the state of the spin. It is 
convenient to choose the unit vectors in such a way that they point into triangles 
of one orientation (say. A). The Hamiltonian of the model with nearest-neighbor 
(nn) exchange Ji can be written in two ways: 

Fi = -Ji^Si.S, = ^5^a.a„ (2.1) 

nn nn 

where we relied on the result • = — 1/2 for nearest neighbors. As usual, a 
ferromagnetic exchange Ji > for spins Sj translates into an antiferromagnetic 
interaction for the Ising variables a^. It is well known that the Ising antiferromagnet 
on kagome remains paramagnetic down to zero temperature (Syozi 1951). It retains 
a finite entropy density in the ground state, 0.502 per spin (Kano & Naya 1953). 

To understand this residual spin degeneracy, we rewrite Hi in terms of magnetic 
charges of triangles 

Hi = ^EQI- (2-2) 

a 

Here we have neglected an irrelevant constant. The magnetic charge is defined as 

0„ = ±^a„ (2.3) 
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with the plus sign for one type of triangles and minus for the other. As discussed 
above, the allowed values of magnetic charges are Qa — ±3 and ±1. Eq. (j2.2p 
shows that triangles with triple charges are energetically disfavored. Consequently, 
as temperature is lowered the magent gradually enters the spin-ice phase consisting 
exclusively of triangles with = ±1, corresponding to the two-in-one-out and 
one-in-two-out ice rules on kagome. The number of the spin-ice microstates grows 
exponentially with the number of spins. Invoking the famous Pauling estimation 
gives an entropy density Sicc — (1/3) ln(9/2) — 0.5014 per spin, which agrees very 
well with the numerical results. 

Wills et al. (2002) added interactions between next- nearest neighbors (nnn): 

= - J2 E • = Y E ^''^^ • (2-4) 

nnn nnn 

Whereas they considered both signs of J2 , we are interested in the antiferromagnetic 
case J2 < 0. MoUer and Moessner (2009) and our group (Chern et al. 2011) added 
long-range dipolar interactions, 

H, = ^l^a.<j, ^— , (2.5) 

where D — /io/i^/(47rr^jj) is the characteristic strength of dipolar coupling, are 
spin locations, f y = (r; — Vj)/\vi — rj|, /i is the magnetic dipolar moment of a spin 
and is the distance between nearest neighbors. For brevity, we shall refer to the 
model of Wills et al. (2002), with the Hamiltonian Hi + H2 and antiferromagnetic 
J2 < 0, as the short-range kagome ice and call the model with dipolar interactions, 
with the Hamiltonian Hi + H^, the long-range kagome ice. 

(b) The ground states 

In both models, interactions between Ising variables ai are antiferromagnetic 
for nearest neighbors and ferromagnetic for next-nearest neighbors. On the basis 
of that, one might expect the same type of magnetic order at low temperatures. 
The ground states for the Ising model on kagome with first and second-neighbor 
interactions are known exactly (Wolf & Schotte 1988; Takagi & Mekata 1993). They 
have an enlarged {^/3 x ^/3) magnetic unit cell with 9 spins. The corresponding 
ground states of the ice model are shown in Fig. [UJa). With the aid of a complex 
order parameter, 

^E(T,exp(iQ-r,), (2.6) 

i 

where Q — (47r/3, 0) and N is the number of spins, we can see that they are similar 
to those of the six-state clock model. The phase of the order parameter M = |M|e"^ 
takes on the six values (f> = mr/3, where n is an integer. Fig. [2lja). 

We were able to show that the long-range model with dipolar interactions has 
the same ground states. To do so, we treated the Hamiltonian Hi+H^ = aiAijaj/2 
as a quadratic form in the Ising variables. We replaced the constraint = ±1 with 
a less stringent normalization condition, aiCTi = N, where N is the number of sites, 
and minimized the quadratic form by finding the lowest (most negative) eigenvalue 
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Figure 2. (a) The ground states of the dipolar spin ice on kagome have a complex order 
parameter M = |M|e"* with (j) — mv/3. (b) A proposed intermediate phase (Takagi & 
Mekata 1993; Wills et al. 2002) has tf) = {n + l/2)n/3. In both cases, the magnetic unit 
cell contains 9 sites. 



of the matrix Aij . Although this procedure minimizes the energy over an enlarged 
space of states, the method gives the right answer if the eigenstate with the lowest 
eigenvalue has ai = ±1. That is indeed the case for the Hamiltonian Hi + H^- The 
corresponding eigenstates are the six ground states of the Ising model on kagome. 

(c) Phase transitions: general considerations 

Previous theoretical studies have revealed that the six-state clock model in two 
dimensions may order in a number of scenarios (Cardy 1980). In all of them, there is 
a high-temperature disordered phase with (M) — and a low-temperature ordered 
phase with (M) ^ 0. 

(A) Two Kosterlitz-Thouless (KT) transitions with an intermediate critical phase 
in which (M) = and spin correlations decay as a power of the inverse 
distance. 

(B) From the paramagnetic phase, the system enters a partially ordered phase via 
an Ising transition. A second transition of the 3-state Potts universality class 
takes the system to the ordered state. The intermediate phase has an Ising 
order parameter (M^) ^ 0, whereas (M) — 0. 

(C) Similar to (B) but with the Ising and Potts transitions exchanged. The in- 
termediate phase has a 3-state Potts order parameter (Af^) ^ 0, whereas 
(M) = 0. 

(D) A discontinuous phase transition between the paramagnetic and fully ordered 
phase. 

Numerical studies (Wolfe & Shotte 1988; Takagi & Mekata 1993) provide com- 
pelling evidence for scenario (A) with two KT transitions in short-range Ising mod- 
els on kagome. Takagi & Mekata mistakenly ascribed a nonzero order parameter 
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(M) to the intermediate phase, which is expected to have quasi-long-range order 
with power-law spin correlations and (M) = 0. 

Wills et al. (2002) also suggested an intermediate ordered phase with six possible 
states shown in Fig. [2jb) . They are equivalent to the intermediate states of Takagi 
& Mekata and have a nonzero order parameter (M) with (j) = (n-|-l/2)7r/3. Such an 
intermediate state is clearly inconsistent with any of the four scenarios mentioned 
above. 

The intermediate critical phase in scenario (A) is most simply understood in 
the continuous version of the clock model, namely the XY model with a sixfold 
anisotropy term ag cos 6^ in the free-energy functional. This term is irrelevant in 
the renormalization-group sense above the lower phase transition, so that the inter- 
mediate phase behaves as if the sixfold anisotropy were absent. In principle, there is 
a possibility that the sixfold anisotropy changes sign: at low temperatures, ag < 
yields the states shown in Fig. [2Ia), whereas at higher temperatures, ag > sta- 
bilizes the states shown in Fig. [2jb). However, if the intermediate phase is indeed 
ordered and fully breaks the Zq symmetry, the transition between that phase and 
the fully symmetric paramagnet should follow one of the above scenarios, which 
requires yet another intermediate phase or a discontinuous phase transition. We see 
no evidence for either possibility. 

Our proposal of the intermediate state with ordered magnetic charges corre- 
sponds to scenario (B). The charge-ordering transition is of the Ising type. It takes 
the system from the paramagnetic phase into a state where the Ising order pa- 
rameter (M^) has a nonzero expectation value. If (M^) = (iMl^e"^"^) > 0, the 
state can be thought of as a mixture of the clock states with (p = 0, 27r/3, and 
— 27r/3, Fig. [21a). At a lower temperature, the system selects one of these states in 
a transition of the 3-state Potts universality class. 

To test which of the scenarios are realized in kagome spin ice, we have performed 
large-scale numerical simulations of both the short and long-range models with the 
Hamiltonians Hi + H2 and Hi + H^- 



3. Numerical simulations 

(a) Kagome ice with short-range interactions 

Here we provide numerical evidence for scenario (A), namely two sucessive KT 
transitions, in the short-range ice model on kagome. A standard Metropolis im- 
portance sampling was used to simulate the behavior of the model on L x L unit 
cells with periodic boundary conditions; the number of spins N — 3L^. In a sin- 
gle Monte Carlo step, each spin in the lattice is updated once with a probability 
p = min(l, exp(— Ai^/fc^T)), where AE is the energy change of flipping a spin. It 
is worth noting that the single-spin update is sufficient for both the spin-ice regime 
and even the intermediate critical phase; the corresponding acceptance rates are 
roughly 30 % and 15 %, respectively. About 10000 initial Monte Carlo sweeps 
were discarded to allow the system to equilibrate. We then retained results from 
10^ ^ 10^ Monte Carlo sweeps for computing averages. Our most complete data 
set was taken with parameters J2 = — Ji/3, for which we describe our analysis in 
detail in the following. 
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Figure 3. Temperature dependence of (a) Specific heat C = {{E^) - (Ef) /NksT^, (b) 
magnetic diarge Q, and (c) magnetic order parameter M for varying lattice sizes obtained 
from Monte Carlo simulations of the short-range spin ice model. The simulation was done 
with J2 = — Ji/3. The shaded area indicates the intermediate critical phase. The transition 
temperatures Td = 0.735Ji and Tc2 = 0.845Ji are determined using the finite-size scaling 
relation Eqs. (|3.4p . 




Figure 4. The log-log plots of order parameter M as a function of system size L at various 
temperatures. The solid curves indicate the linear behavior which corresponds to a pow- 
er-law dependence, M oc L-'"'^ in the intermediate critical phase. The inset shows the 
extracted critical exponent r; as a function of temperature in the critical regime. The two 
data points represented by symbol □ are obtained using the finite-size scaling Eqs. (|3.4|l . 
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Figure 5. Distribution of instantaneous magnetic order parameter M in the complex 
plane (ReM, ImM) for different temperatures. 

Fig. [3] shows the temperature variation of specific heat, staggered charge order, 
and magnetic order parameters for different lattice sizes. While the appearance of 
two peaks in the specific-heat curve clearly indicates two phase transitions, the 
locations and amplitudes of the peaks vary only slightly with lattice size, a feature 
characteristic of KT transitions. The intermediate critical phase determined by 
finite-size scaling to be discussed below is marked by the shaded region in the 
figures. The charge order parameter is defined as difference of magnetic charges on 
the two different types of triangles: 

Q = ^E(-1)"^?- (3.1) 

where (—1)" = ±1 for up and down triangles, respectively, and 2L^ is the num- 
ber of triangles in the lattice. The charge order is proportional to the Ising order 
parameter, (Q) c>c (M^). As can be seen from Fig. [3Jb), the staggered charge or- 
der decreases rapidly with increasing L above the low-T transition, indicating that 
magnetic charges remain disordered in the intermediate phase. 

On the other hand, the magnetic order parameter shows quite different finite- 
size behavior across the three different regimes. Fig. [3jc). Above the high-T peak 
the rapid decrease of M with increasing L is consistent with a magnetically disor- 
dered phase, whereas a negligible finite-size effect implies an ordered state below 
the second transition at Td. In the intermediate regime the order parameter M falls 
off rather slowly with increasing system size. These observations are summarized 
in Fig. |4] which shows magnetic order parameter M as a function of L at different 
temperatures. The magnetic order parameter extrapolates to zero at high temper- 
atures, whereas it levels out to a constant at low T. In the intermediate regime, the 
linear behavior in the double-logarithmic plot indicates a power-law dependence 

McxL-'''/^ (3.2) 

which is consistent with a divergent correlation length in a critical phase. The 
extracted value of critical exponent 77 as a function of temperature is shown in the 
inset of Fig. H] 

As discussed in Sec. [HeI), the sixfold anisotropy gq cos Gcj) in the Landau expan- 
sion of the free-energy functional is irrelevant in the intermediate critical phase. To 
confirm this, we show in Fig. [S] the distribution of instantaneous order parameter 
M = |M|e"^ in the complex (ReM, ImM) plane at various temperatures. As can be 
seen in Fig.jSKc) which corresponds to a temperature in the intermediate regime, the 
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Figure 6. (a) Finite-size scaling of the order parameter M for the lower-temperature tran- 
sition at Tel. (b) Finite-size scaling of susceptibihty xm = {{M'^) — {M)^)/NkBT for the 
higher-temperature transition at Tc2- The following values of the parameters are used: 
a = 1.22, b = 0.0615, c = 1.746, Td = 0.735 Ji, and Tc2 = 0.845 Ji. 



order parameter has a finite amplitude \M\ ^ due to finite system size (L = 180). 
More importantly, a continous 0(2) symmetry emerges in this critical phase. In the 
ordered phase below Td, the anisotropy term reasserts itself and restores the Zq 
symmetry as demonstrated in Figs.lSJa) and (b). 

To further corroborate the proposed scenario of two KT transitions, we resort 
to the method of finite-size scaling (Challa & Landau 1986). For a KT transition, 
the correlation length near the critical temperature diverges as (Kostcrlitz 1974) 

^ cx exp(at"i/2)^ (3 3) 

where a is a non- universal constant and t ~ \T — Tc\/Tc is the reduced temperature. 
The order parameter and its susceptibility exhibit power law behavior, M oc ^"''/^ 
and Xm oc in a KT transition. For finite systems the singular part of the free 

energy is expected to depend only on ^/L. The order parameter and susceptibility 
are assumed to have the functional forms 

M = L-^/^M{^/L), XM = i'-" Xi^L), (3.4) 

where A4 and X are unknown universal functions. At the critical point, substituting 
^ — oo into Eq. p. 41) gives the power law relation p.2p . For an extended regime 
consisting of a line of critical points, the power law M ~ L"''/^ should hold over 
the entire temperature range from Td to Tc2, which is indeed observed in our 
simulation. Fig. HI 

The finite-size scaling relations p.4p can also be used to determine the upper 
and lower critical temperatures of the intermediate KT phase. For example, with 
an appropriately chosen set of parameters a, b, and Td, the plot of ML'' versus 
exp{a/^/Tci — T) should collapse on a universal curve for different system sizes. 
The same should hold for xm-^ versus exp{a/^/T — Tc2)- The two constants, b 
and c, are related to the critical exponent: b — 7//2 and c — 2 — 7y. As shown in Fig. El 
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Figure 7. Temperature dependence of (a) Specific lieat C = {{E^) - {Ef)/NkBT^, (b) 
charge Q, and (c) magnetic M order parameters for varying system sizes. The simulation 
was done with D = 2 Ji . 

excellent data collapse was obtained using the following set of parameters: a = 1.22, 
b = 0.0615, c = 1.746, Tci = 0.735Ji, and Tc2 = 0.845Ji. From this, we estimated 
the critical exponent at the two transition temperatures: r]{Tci) = 2b = 0.123 and 
v{Tc2) = 2 — c = 0.255. These values agree resonably well with those extracted by 
linear-fitting of the power-law relation (I3.2p within the critical phase, see the inset 
of Fig. m Both values are close to the theoretical predictions 1/9 and 1/4 for the 
six-state clock model (Jose et al. 1977), and the slight overestimation could be due 
to the finite-size effects on KT transitions. 

(&) Kagome ice with dipolar interactions 

Although the single-spin Metropolis algorithm is quite efficient for simulating 
short-range kagome ice, it suffers from a dynamical freezing in the charge-ordered 
phase. A similar problem arises in the low-temperature simulation of pyrochlore spin 
ice (Melko et al. 2001) in which the single-spin flip violates the ice rules, leading to a 
large energy cost and low acceptance rate of the updates. To overcome this problem 
we added the nonlocal loop moves first introduced by Barkema & Newman (1998) 
for square ice models. In a nonlocal update, a loop is first formed by randomly 
tracing a path through triangles satisfying the ice rules, alternating between spins 
pointing in and spins pointing out of the triangles. The process is completed when 
the path closes upon the starting spin or encounters any spin already included in 
the loop. Flipping all the spins in such a loop leaves the magnetic charges {Qa} 
intact and conserves the nearest-neighbor energy Hi. The loop move results in a 
small gain or loss of the dipolar energy A£^d, and the update is accepted with a 
probability p = min(l, exp(— Ai^d/fesT"))- 

We employed a combination of single-spin flips and loop moves to simulate 
the long-range ice model Hi + H^. With periodic boundary conditions, dipolar 
interactions were summed over periodic copies up to a distance of 500L. Most of 
the results presented here are obtained with a dipolar coupling D = 2Ji. Fig. [7] 
shows the temperature dependence of specific heat, charge and magnetic order 
parameters for various system sizes. Similar to the case of short-range kagome ice. 
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Figure 8. Finite-size scaling of (a) Binder's fourth-order cumulant _B4q — 1 — (Q'')/3{Q^)^, 
(b) specific heat, (c) staggered charge order parameter, and (d) charge susceptibil- 
ity XQ = ((Q^) ~ {Q)^)/NkBT for the high-temperature transition. The inset in (a) 
shows Binder's cumulant of different lattice sizes crossing at the critical temperature 
Tc2 ~ 0.29-D. The critical exponents of two-dimensional Ising model a = 0, /3 = 1/8, 
7 = 7/4, and u — 1 are used. 



two peaks can be seen in the specific-heat curves, indicating two phase transitions. 
While the high-temperature peak becomes sharper with increasing L, the low-T 
transition barely shows any finite-size dependence. The behaviors of charge and 
magnetic order parameters shown in Figs. [Tj^b-c) are consistent with scenario (B) 
discussed in the Introduction. We discuss both transitions below. 

The staggered charge order corresponds to the partially ordered phase in sce- 
nario (B) which is characterized by an Ising order parameter (Q) (x (M^) . The 
charge-ordering transition is thus expected to be in the universality class of the Ising 
model. To verify this conjecture, we performed finite-size scaling analysis and found 
excellent data collapse. Fig. [8j using the critical exponents of the two-dimensional 
Ising universality class. The critical temperature Tc2 = 0.29D is determined from 
the crossing of the Binder's fourth-order cumulant for different lattice sizes. 

The origin of an intermediate phase with ordered charges can be understood 
by expressing the energy of the dipolar ice in terms of magnetic charges. This is 



Article submitted to Royal Society 



Kagome spin ice 



13 



C 

0.14 
0.12 
0.1 
0.08 
0.06 



















/ 








L=9 — ' — 
12 

15 --«- ■ 

24 a 

36 - - - 





0.12 



0.16 

TjD 



0.12 



0.08 



0.2 




0.04 



Figure 9. Specific iieat as a function of temperature for (a) dipolar ice model on kagome, 
and (b) dimer model with next-nearest-neiglibor attraction on tlie liexagonal lattice. The 
inset in (b) shows the temperature variation of the order parameter characterizing the 
\/3 X \/3 dimer order. 



achieved through the so-cahed dumbbeh approximations (Castelnovo et al. 2008) 
in which the dipoles are stretched into bar magnets of length a — (2/\/3)rnn such 
that their poles meet at the centers of triangles: 

Here = /i/a is the magnetic charge of the dumbbell, ^ is the moment of the 
spin. The self-energy vq = J/2 + (11 + 3%/3)D/8 for ka gome ice. It can be shown 
that the dipolar energy Hi + is well approximated by Eq. p. 51) with corrections 
which vanish with distance at least as fast as for each dipole pair. 

The dumbbell approximation also illustrates an important difference bwteen 
pyrochlore and kagome ices. Since the allowed charges are even on a tetrahedron 
and odd on a triangle, minimization of the self-energy term results in Qa = on the 
pyrochlore lattice and ±1 on kagome. The conditions of minimal allowed charges on 
simplexes correspond to the ice rules on the respective lattices. For pyrochlore ice, 
the vanishing of magnetic charge also makes the Coulomb interaction, the second 
term in Eq. p.Sp . strictly zero. This observation explains why all states satisfying 
the ice rules are esstentially degenerate in pyrochlore spin ice over a wide range of 
temperatures (Gingras & den Hertog 2001). The degeneracy is lifted by the residual 
interactions neglected in the Coulomb approximation and a long-range magnetic 
order with wavevector Q — (0, 0, 27r) sets in at a lower temperature compared to 
the dipolar energy scale (Melko et al. 2001). 

The situation on kagome is quite different as nonzero charges generate a mag- 
netic field. This results in substantial energy differences between states obeying 
the ice rule Qa = ±1- Consequently, interactions between uncompensated charges 
induce an Ising transition into a state with staggered charge order which minimizes 
the Coulomb energy. This partially ordered phase is closely related to the ice states 
of pyrochlore spin ice in a (111) magnetic field (Udagawa et al. 2002; Moessner & 
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Figure 10. Top panel: spin configurations. Bottom: respective dimer configurations, (a) 
Reference state witii ciiarges Qa = ±3 on the two sublattices. (b) A charge-ordered state 
with Qa = ±1. (c-e) two charge defects {AQa — ±2) are created and pulled apart. Red 
arrows denote the minority spins. 



Sonhdi 2003). Triangles with positive (negative) charges on kagome correspond to 
three-in-one-out (three-out-one-in) tetrahedra of the pyrochlore lattice. 

Spins remain disordered in the charge-ordered ice phase as evidenced by the 
loop moves discussed at the beginning of this section. Such nonlocal loop updates 
flip spins on a closed path while maintaining the charge configuration. To quantify 
this residual degeneracy, we note that each triangle in a state with perfect charge 
order has two majority spins pointing into (or out of) the triangle and a minority 
spin pointing the other way. Such states can be mapped to dimer coverings on the 
honeycomb lattice by identifying the minority spins with the dimers, Fig. llOf b). 
The number of dimer-coverings on honeycomb grows exponetially with the lattice 
size, giving rise to a residual entropy density S ~ 0.108 per spin (Udagawa et al. 
2002). 

The remaining entropy of the charge-ordered phase is completely removed by 
a magnetic phase transition corresponding to the low-T peak of the specific-heat 
curve, Fig. EKa). The magnetic order is shown in Fig. [Ud) and is again charac- 
terized by order parameter M; it has an enlarged a/S x \/3 unit cell containing 9 
spins. Similar to pyrochlore spin ice, the magnetic ordering is induced by residual 
interactions beyond the dumbbell approximation Eq. (j3.5p . 

As the Zq symmetry of the spin-ice Hamiltonian is reduced to a threefold sym- 
metry in the charge-ordered phase, the magnetic transition is expected to be in 
the universality class of 3-state Potts model. However, Monte Carlo simulations on 
systems up to L = 36 fail to turn up any signature of the Potts criticality. Instead, 
the lack of a singularity in the specific heat, Fig.j^Ja), seems to be consistent with 
a KT transition. 

(c) Dimers on a honeycomb lattice 

To shed light on this puzzling result we consider a similar transition in the hon- 
eycomb dimer model. As discussed above, the charge-ordered states can be uniquely 
mapped to dimer coverings on honeycomb. To induce the corresponding ^73 x ^73 
ordering of dimers, we introduce a next-nearest-neighbor attractive interaction be- 
tween the dimers. Note that nearest-neighbor dimers are precluded by the hardcore 
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constraints that each lattice site is attached to exactly one dimer. The partition 
function of the modified dimer model reads 

„ >^ fvMcy\ 

where the sum is over all dimer coverings of the hexagonal lattice and N2{C) counts 
the number of next- nearest- neighbor pairs in a given covering C. An attractive 
interaction between next-nearest-neighbor dimers corresponds to V2 > 0. 

We employed the direct-loop Monte Carlo algorithm (Sandvik & Moessner 2006) 
to simulate the dimer model p.6p . The nonlocal update is performed by initially 
breaking up an arbitrary dimer into a pair of monomers and then moving one 
monomer across the lattice by flipping a sequence of dimers along a path. The pro- 
cess is completed and a new dimer configuration is created once the two monomers 
meet and recombine with each other. Since detailed balance is maintained at each 
step of the monomer's movement, large number of dimers can be efficiently updated 
without rejection. This method significantly reduces the autocorrelation time in the 
Monte Carlo process and allows us to simulate large lattices up to L = 150. Fig.lHI^b) 
shows the specific heat versus temperature for various system sizes; also shown in 
the inset is the temperature dependence of the order parameter characterizing the 
1/3 X -x/S dimer order. The rather weak finite-size dependence and a lack of sin- 
gularity in specific heat again shows that the dimer ordering is characterized by a 
KT transition similar to the case of dimer model with aligning interactions on the 
square lattice (Alet et al. 2006). 

The appearance of a KT transition indicates the critical nature of the disor- 
dered hardcore dimer phase, or equivalently the charge-ordered ice phase. Indeed, 
by further mapping the dimer covering to a 'height' field (Blote & Hilhorst 1982), 
the dimer model is described by a sine-Gordon model in the coarse-grained approx- 
imation. It is well known that the height model undergoes a KT transition into a 
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temperature Tm ~ 1.226i;2. The critical exponents of two-dimensional 3-state Potts 
model a = 1/3, /3 — 1/9, 7 — 13/9, and 1/ = 5/6 are used. 



rough phase at high temperatures (Kardar 2007). As the confining cosine potential 
term is irrelevant at high-T, the rough phase is described by a Gaussian field theory 
with critical correlation functions in two dimensions. 

The existence of such a critical phase relies on the absence of charge defects. At 
finite temperatures, however, thermally excited defects spoil the mapping to hard- 
core dimer coverings. There are two types of defects in the charge-ordered phase: 
triangles with ±3 charges (ice defects) and triangles that violate the charge order 
(charge defects). In the dimer model, the triply charged sites become monomers, 
whereas defects in charge order become sites with two dimers. Fig. [TO] We ig- 
nore the triply charged defects in the intermediate phase on account of their high 
energy cost and focus on the charge defects. To that end, we modified the honey- 
comb model (j3.6l) by allowing dimer pairs with fugacity z = exp(— e/fc^T). The 
hardcore-dimer model is recovered in the limit e — 00. Fig. [TT] shows the tempera- 
ture variation of the a/S x a/3 order parameter, susceptibility and density of dimer 
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pairs for e = 2, 10 and oo in units of V2- The model with larger fugacity for dimer 
pairs shows a dramatically different behavior from that of hardcore dimer covering. 

For e < +00, a finite density of dimer pairs sets a bound on the dimer corre- 
lation length, thus altering the criticality of magnetic ordering to the universality 
class of 3-state Potts model. To confirm this, we performed a finite-size scaling 
study on the dimer model with e — 2v2 ■ As can be seen from Fig. [TTJ one obtains 
excellent data collapse with the critical exponents of the two-dimensional 3-state 
Potts model. However, when the average separation between charge defects exceeds 
the lattice size, we are back to hardcore dimers with power-law spatial correlations 
for all distances. This explains the KT-like behavior observed in spin ice at small 
lattice sizes. The critical behavior characteristic of the 3-state Potts universality 
only reveals itself for sufficiently large systems (Otsuka 2011). 

4. Discussion 

Magnetic charges in spin ice are an example of an emergent phenomenon. Although 
the fundamental degrees of freedom in these materials are magnetic dipoles, the 
behavior of the system at low energies is often more conveniently expressed in the 
language of magnetic charges. In spin ice on kagome, magnetic charge of a simplex 
can take on odd values ±1 and ±3 in natural units. Even though magnctostatic 
interactions tend to minimize magnetic charge, low-energy spin-ice states have non- 
vanishing magnetic charges on simplexes. Such a magnet may have, in addition 
to the paramagnetic and fully ordered states, a distinct intermediate phase with 
ordered magnetic charges and disordered spins. This phase possesses considerable 
residual entropy (0.108 per spin) because a given pattern of magnetic charges can 
be realized by many configurations of magnetic dipoles. 

Our numerical simulations provide evidence for an intermediate phase with stag- 
gered magnetic charges in spin ice with dipolar interactions on kagome. As the 
system is cooled down from the paramagnetic state, it first gradually enters the 
spin-ice regime in which magnetic charges of triangles are restricted to the smallest 
(in magnitude) values of ±1. A phase transition of the Ising universality class takes 
the magnet into the intermediate phase with staggered magnetic charges. 

At a lower temperature, the system enters a magnetically ordered phase with 
six ground states shown in Figure [H^a) . Because only three of them are accessible 
from each of the two states of the charge-ordered phase, this transition is expected 
to be in the universality class of the 3-state Potts ferromagnet. However, observ- 
ing the corresponding critical behavior proved challenging because the transition 
takes place in the background of nearly-perfect magnetic charge order. When the 
charge order has no defects at all, magnetic configurations can be mapped onto 
states of hard-core dimers with attraction. This mapping establishes that spins in 
the background of perfectly ordered magnetic charges are in a quasi-ordered phase 
with algebraic spatial correlations. The transition to the fully ordered phase is then 
of the Kosterlitz-Thouless type. Thermal fluctuations generate defects in charge or- 
der, thus violating the hard-core constraint for dimers. The spin correlation length 
is then set by the typical distance between isolated charge defects. As Figure [7Kb) 
shows, the charge order parameter quickly approaches saturation upon cooling be- 
low Tc2 = 0.291?. At T = 0.24Z3, i.e., considerably above the onset of magnetic order 
at Tci, the concentration of charge defects is about 0.01 per triangle. Furthermore, 

Article submitted to Royal Society 



18 



G.-W. Chem and 0. Tchemyshyov 



most of these defects come in pairs with zero net charge and minimal separation (re- 
sulting from single spin flips), so they just renormalize the spin correlation function 
without shortening the correlation length. Not surprisingly, the intermediate phase 
exhibits power-law spin correlations up to very large distances and the transition 
to the magnetically ordered phase appears to be of the Kosterlitz-Thouless type. 
Observing the true critical behavior of the 3-statc Potts universality class would 
require extremely large system sizes inaccessible to us. 
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